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Abstract 

A very fast heuristic iterative method of projection on simphcial cones 
is presented. It consists in solving two hnear systems at each step of the 
iteration. The extensive experiments indicate that the method furnishes 
the exact solution in more then 99.7 percent of the cases. The average 
number of steps is 5.67 (we have not found any examples which required 
more than 13 steps) and the relative number of steps with respect to 
the dimension decreases dramatically. Roughly speaking, for high enough 
dimensions the absolute number of steps is independent of the dimension. 

1 Introduction 

Projection on polyhedral cones is one of the important problems applied op- 
timization is confronted with. Many applied numerical optimization methods 
uses projection on polyhedral cones as the main tool. 

In most of them, projection is part of an iterative process which involve its 
repeated application (see e. g. problems of image reconstruction [1 , nonlinear 
complementarity [3J [H], etc.). Hence, it is important to get fast projection 
algorithms. 
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The main streems of the current methods in use rely on the classical von 
Neumann algorithm (see e.g. the Dykstra algorithm O [U [10] )j but they are 
rather expensive for a numerical handling (see the numerical results in [7] and 
the remark preceding section 6.3 in [5]). 

Finite methods of projections are of combinatorial nature which reduces 
their applicability to low dimensional ambient spaces. 

Recently we have given a simple projection method exposed in note jB] for 
projecting on so called isotone projection cones. Isotone projection cones are 
special simplicial cones, and due to their good properties we can project on 
them in n steps, where n is the dimension of the ambient space. In the first 
part of that note we have explained our approach by considering the problem of 
projection on simplicial cones by giving an exact method based on duality. This 
method has combinatorial character and therefore it is inefficient. More recently 
we observed that a heuristic method based on the same ideas gives surprisingly 
good results. This note describes the theoretical foundation of the heuristic 
method and draws conclusions based on millions of numerical experiments. 

Projection on polyhedral cones is a problem of high impact on scientific 



2 The simplicial cone and its polar 

Let M" be an n-dimensional Euclidean space endowed with a Cartesian reference 
system. We assume that each point of M" is a column vector. 

We shall use the term cone in the sense of closed convex cone. That is, the 
nonempty closed subset K C K" in our terminology is a cone, if K + K (Z 
and tK C K whenever i e M, t>0. 

Let m < n and ei, . . . , e„ be m elements in M". Denote 

cone{ei, . . . , em} — {A^ei + • • • + A™e,„ : A* > 0, i = 1, . . . , m}, 

the cone engendered by ei, . . . , e„i- Then, 

cone{ei, . . . , e„.} = {Ev : v € M™}, (1) 

where E = (ei, . . . , e„j) is the matrix with columns ei, . . . , Bm and is the 
non- negative orthant in M™. 

Suppose that ei, . . . , e„ e M" are linearly independent elements. Then, the 
cone 

K = cone{ei,...,e„} , , 

= {A^ei + • • • + A"e„ : A* > 0, i = 1, . . . , n} = {Ev : v e M™}, 

with E the matrix from ([T]) for m = n, is called simplicial cone. Denote N — 
{1,2,. 

The polar of K is the set 

K° = {x e M" : x^y < 0, Vy G K}. (3) 

K* = —K° is called the dual of K. K is called subdual, if if C K* . This is 
equivalent to the condition e^ej. > 0, £, k G N. 

^see the popularity of the Wikimization page Projection on Polyhedral Cone at 
http: / / www.convexoptimization.com / wikimization/index.php / SpeciahPopularpages. 



2 



Lemma 1 The polar of the simplicial cone can be represented in the form 
K° = {^i^ui + . . . ^"u„ : > 0, i = 1, . . . , n}, (4) 
where Ui(i = 1, . . . , n) is a solution of the system 

ejui = 0, j = 1, . . . ,n, j ^i, 
eju, = -1 

is normal to the hyperplane spanjei, . . . , e^-i, e,;+i, . . . , e„} in the opposite 
direction to the halfspace that contains Gi). Thus, 

K° = {Ux : X e R^} 

with M" = {x = {xi, . . . , : Xi > 0, i — 1 . . .n} and 

U = -(E-i)^. (5) 

For simplicity we shall call U the polar matrix of E. The columns of U are 
{Ui: j = 1, . . . ,n}. 

Proof. Let y = X]j=iM"''^i ^ ~ SiLi o^*^* ^'^^ non-negative real 
numbers a' and /i^ . The inner produce of y and z is non-positive because 

n n n 

y^z = J2T. «V'e7u, = - ^ aV^ < 

i—l j — 1 i—1 

But y is an arbitrary element of the right hand side of Q and z is an arbitrary 
element of K, thus we can conclude that the righ hand side of Q is a subset of 
K°. 

The vectors Ui, . . . , u„ are linearly independent. (This can be verified by 
assuming the contrary, and by multiplying the subsequent relation by to get 
a contradiction.) Hence for y G K° we have the representation 

2/ = ^'ui + ...+/3"u„. 

By (§, 

ejy = 

so f^^ > which prove that y is an element of the right hand side of Q . Thus 
we can conclude that K° is a subset of the right hand side of Q . □ 

The formula ^ of Lemma [l] is equivalent to the formula (380) of T. 

Corollary 1 For each subset I of indices in N, the vectors ei,i G I,Uj,j G /"^ 
( where I'^ the complement of I with respect to N ) are linearly independent. 

Proof. Assume that 
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for some reals a' and (3^ . By the mutual orthogonality of the vectors , i G / 
and Uj,j S /'^ it follows, by multiplication of the relation (joj) with X^ie/'^'^^ 
and respectively with J2jei'' /^"''^J' ^^^^ 

J2 a'e, = 

and 

Hence, a' = /J-' = must hold. □ 

The cone Kq C K is called a face of /-iT if from x e i^'o, y ^ K and x — y e iiT 
it follows that y e iiTo. The face Ko is called a proper face of A', if Kq ^ K. 

Lemma 2 If K is the cone ^ and K° is the cone Q), then for every subset 
of indices I — {ii, i^} C N the set 

Fi = cone{ei : i G /} = {x G if : x^Uj = : i T} (7) 

fzwit/i Fi = {0} if 1 ^ is a face of K . If ih 7^ ii whenever h ^ I, then Fj is 
for k > a nonempty set in M" of dimension k. (In the sense that Fj spans a 

subspace o/M" of dimension k.) 

Every face of K is equal to Fj for some I d N. If I ^ N then Fj is a proper 
face. 

Proof. 

The relation in ^ follows from the definition of the vectors Uj in Lemma 
[l] while the assertion on the dimension of Fj is obvious. 
Suppose that x e F/ and y G if with y < x. 

Then (x - y)^Uj = -y^Uj < 0, Vj e r and y^u^ < 0, Vj € N, because 
y G K. Thus y^Uj = 0, Vj E F, hence y G Fj, showing that Fj is a face. 

Suppose that x G i^ for F an arbitrary proper face of K. Since x G if , by 
the definition of the vectors Uj, x^Uj < for j G N. 

If x^Uj < 0, Vj G N, then there exists a positive scalar t with (x— ty)^Uj < 
0, Vj G N. Hence, x — ty G if and thus ty < x. But then ty G F and since F 
is a cone, y G F. This means that K C F, that is, F cannot be a proper face. 

We have to show that F has a representation like ([t]). By the above reason- 
ing, for each x G F there exist some index i G N with with x^u^ = 0. 

If F = {0} we have the representation ([t]) with i = 0. 

If F 7^ {0}, take x in the relative interior of F and let I be the complement 
in N of the maximal set of indices j with x^u^ — 0. (i must be a nonempty, 
proper subset of N since x 7^ 0.) 

Take y G F arbitrarily. By the definition of x, x— iy G F for some sufficiently 
small t > 0. Hence, 

(x-ty)^u, <0, yieN. (8) 

By y G F C if we also have y^u^ < 0, Vi G N. If y^u^ < for some j G I'^, 
then ^ would imply 

x^Uj < ty^Uj < 0, 
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which is a contradiction. Hence, we must have y Uj = 0, Vj G J'^; and accord- 
ingly 

F C {z e X : z^u^ = 0, Vj e r}. (9) 

Suppose that y ^ K and y^u^ = 0, Vj G I"^. From definition we have 
x^Ui < for each i e /, whereby for a sufficiently large i > 0, 

(tx-y)^u, < 0, Vi e A^. 

Hence, tx— y is in the polar of K° , which by Farkas' lemma is K (This follows in 
fact, in our case, also by the symmetry of the vectors and \ij in the formulae 
of Lemma [l]) Thus < y < tx, whereby < {l/t)y < X. Since is a face 
of K, we have (l/i)y S F and since it is also a cone, y G F. This proves the 
converse of the inclusion in ^ and completes the proof. 

□ 

Thus a maximal proper face of K is of the form 

Ki^ = cone{ei : i e N \ {iq}} = cone{ej : i e N, eju,^ = 0}, 

hence it is also called the face of K orthogonal to u^^ . Similarly, we have a 
maximal proper face of K° orthogonal to some ej„ . 

An equivalent result to the one presented in Lemma [2] is given by the Cone 
Table 1 on page 179 of [T]. 

Thus a maximal proper face of K is of the form 

Ki^ — coneje; : i £ N \ {io}} — conejei : i £ N, ejui^ ~ 0}, 

hence it is also called the face of K orthogonal to u^^ . Similarly, we have a 
maximal proper face of K° orthogonal to some Bj^^ . 

Let F = conejei : i € /} and F^ = conejuj : j G I'^}. Then, from the 
above results it follows that 

F = {xe K : x^Uj = 0, j e F} 

and 

F^={yeK°: y^e, - 0, z G /}. 

The faces F C K and F-^ C K° of the above form are called a pair of 
orthogonal faces where F-^ is called the orthogonal face of F and F is called the 
orthogonal face of F-^ . 

3 Finite method of projection on a simplicial 
cone 

For an arbitrary u S M" denote ||u|| = V u^u. Let K G M" be an arbitrary 
cone and K° its polar, and C C M" an arbitrary closed convex set. Recall that 
the projection mapping Pc : H ^ H on C is well defined by Pcx € C and 

||x - Pcx|| = min{||x - y|| : y G C}. 

Then, Moreau's decomposition theorem asserts: 
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Theorem 1 (Moreau, l^) For:x., y, z G M" the following statements are equiv- 
alent: 

(i) z ^ x + y,x e K,y e K° and x^y = 0. 

(ii) X = P^z and y = Px°z. 

Suppose now, that K is a shiipHcial cone in M". We shall use the represen- 
tation ([2| for K and the representation Q for K° . Hence, 



e7u, 



-S),i,j = 1, 



where the Kronecker symbol. As a direct implication of Moreau's decompo- 
sition theorem and the constructions in the preceding section we have: 

Theorem 2 Let x G M" . For each subset of indices I C N , x can be represented 
in the form 

x = 5]a*e,+ ^/3^u, (10) 

iei jei" 

with /"^ the complement of I with respect to N , and with a* and (3^ real numbers. 
Among the subsets I of indices, there exists exactly one (the cases / = and 
I = N are not excluded) with the property that for the coefficients in (10) one 
has /?■' > 0,j G /'^ and a* > 0,i G /. For this representation it holds that 

PA'X==^a*e„ a'>0, (11) 



lei 



and 



> 0. 



(12) 



Proof. The first assertion is the consequence of Corollary [T] 
The projections P/fX and PifoX as elements of K and K° , respectively can 
be represented as 



P^.x = ^a*e„ a'>0 



4=1 



and 



(13) 



(14) 



To prove existence, let I'^ = {j E N : (5^ > 0} and let / be the complement of /"^ 
in the set N of indices. For an arbitrary element z G M", denote P]^ z ~ (P^f z)^. 
If > would hold in (14), for some j G I'^, then by Lemma [l] it would follow 
that Pjx-Pifox < 0, which contradicts the theorem of Moreau. Hence, ( 13 ) can 



be written in the form (111 and ( 14 1 can be written in the form ( |12[ ). Therefore, 
Theorem [T] implies 

X = Pa'X + Pifox = ^ a'ei + ^ 



where a* > 0, Vi G / and (3^ > 0, Vj G 
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To prove uniqueness, suppose that in the representation ( 10 ) of x we have 
a* > 0, = for i G I and > 0, = for j S I'^, where / is a subset 
of N, and I'^ is the complement of / in iV (the cases I — % and I — N are not 
excluded). Then representations (11) and ( [l2| follow from Theorem [l] by using 
the mutual orthogonality of the vectors e,, is/ and Uj, j G 7*^. From (11) 
and the uniqueness of the projection Pi^x it follows that / is unique. 

□ 

From this theorem it follows that a given simplicial cone K C M" determines 
a partition of the space K" in 2" cones in the sense that 



(J cone{ej,Uj : i£ I, j e r} 



ICN 

and for two different sets / of indices the respective cones do not contain common 
interior points. The cones in the above union are exactly the sums of orthogonal 
faces. 

This theorem suggests the following algorithm for finding the projection 

Step 1. For the subset I C N we solve the following linear system in a* 

x^e£ ^ a^ejee, I e /. (15) 
iei 

Step 2. Then, we select from the family of all subsets in N the subfamily A 
of subsets / for which the system possesses non-negative solutions. 
Step 3. For each / e A we solve the linear system in 

x^Ufc = /^'ujufc, k e r. (16) 

By Theorem [2] among these systems there exists exactly one with non-negative 
solutions. By this theorem, for corresponding / and for the solution of the 
system ( |15[ ), we must have 

P^x = ^a'ei. 

is/ 

This algorit hm requires that we solve 2" linear syste ms of at most n equa- 



tions in Step 1 (15) and another 2l^l systems in Step 2 (16). (Observe that all 
these systems are given by Gram matrices, hence they have unique solutions.) 
Perhaps this great number of systems can be substantially reduced, but it still 
remains considerable. 

Remark 1 If K is subdual; that is, ifejei > 0, fc, ^ G N, the above algorithm 
can be reduced as follows: By supposing that we have got the representation H 



of X with non-negative coefficients, we multiply both sides of {10) by an arbitrary 
ej . Ifx^ei < then I cannot be in I, otherwise the relations Ujei = 0, j E I'^ 
and ejei > 0, i E I would furnish a contradiction. Thus, we have to look for 



the set I of indices (for which we have to solve the system (15)) among the 
subfamilies of {i E N : x^e^ > 0}. (Arguments like this can be used, as it was 
done e. g. in ^ for the Dykstra algorithm, to eliminate some hyperplanes while 
comkputing successive approximations of the solution.) 
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Obviously, the proposed method is inefficient. It was presented by A. B. 
Nemeth and S. Z. Nemeth in 18] as a preparatory matherial for an efficient 
algorithm for so called isotone projection cones only. For isotone projection 
cones we can obtain the projection of a point in at most n steps, where n is the 
dimension of the space. Isotone projection cones are special simplicial cones. 
Even if there are important isotone projection cones in applications, they are 
rather particular in the family of .simplicial cones. 

4 Heuristic method for projection onto a sim- 
plicial cone 

Regardless the inconveniences of the above presented exact method, which follow 
from its combinatorial character, it suggests an interesting heuristic algorithm. 
To explain its intuitive background we consider again the simplicial cone 

K = cone{ei, . . . ,e„} 

and its polar 

K° =cone{ui,...,u„} 

given by Lemma [l] 

Take an arbitrary x €E M". We are seeking the projection Pkx. 

If e^x < 0, V i e iV, then x G K° = kerP^f , hence P^x = 0. 

If ujx < 0, y j ^ N, then x G K, and hence P_r-x = x. 

We can assume that x ^ K U K° . Hence, x projects in a proper face of K 
and in a proper face of K° . 

Take an arbitrary family I C N oi indices. Then, the vectors 

e^, : ie I, j e r 
ent gender by Corollary [T] a reference system in M". Then, 

x = ^a*ei+ ^/S^Uj (17) 

iei jei" 

with some a', f3^ G R. (As far as the family / C iV of indices is given, we can 
determine the coefhcients a* and (3^ , according to Theorem [2] by solving the 



systems (15) and (16).) 

If we have a* > 0, > : i E I, j E then from Theorem [2] we obtain 



a Bi 

iei 



and 



jei'= 

In this case x is projected onto face F — cone{e,i : i E 1} ortogonally along 
the subspace engendered by the elements {uj : j E I'^}, roughly speaking, along 
the orthogonal face F-^ of F. 

Suppose that < for some j E I'^. Then, considering the reference system 
entgendered by e^, Uj, i E I, j E I'^, x lies in its orthant with negative j*^ 
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coordinate, that is in the direction of the vector —uj. By construction, and 
Uj form an obtuse angle. Hence the angle of and —Uj is an acute one. Thus 
there is a real chance that in a new reference system in which Bj replaces Uj, 
the coordinate of x with respect to has the same sign as its coordinate with 
respect to — u^, that is positive (or at least non- negative) . 

If we have a' < for some i € /, then by similar reasoning it seems to be 
advantageous to replace with u^, and so on. 

Thus, we arrive to the following step in our algorithm: 

Substitute Uj with if f3^ < and substitute with if a' < and solve 
the systems ( 15 ) and ( 16 ) for the new configuration of indices /. We shall call 
this step an iteration of the heuristic algorithm. 

Then, repeat the procedure for the new configuration of / and so on, until 
we obtain a representation (17) of x with all the coefficients non-negative. 



5 Experimental results 

The heuristic algorithm was programmed in Scilab, an open source platform for 
numerical computationj^ Experiments were performed on numerical examples 
for 2, 3, 5, 10, 15, 20, 25, 30, 50, 75, 100, 200, 300, 500 dimensional cones. The 
algorithm was performed on 100000 random examples for each of the problem 
sizes 2, . . . , 100. Statistical analysis on a subset of 10000 examples from the 
set of 100000 examples for size 100 indicates no significant difference in overall 
results and performance, therefore we subsequenlty reduced the number of ex- 
periments on larger problem sizes. 10000 random examples were used for sizes 
200 and 300 and 1000 examples for size 500, as the time needed by the algorithm 
increases with size. Table [T] shows the experimental results. For each problem 
size, the averages of all runs are shown, together with a confidence interval at 
confidence level 95% where appropriatej^ The Changes column indicates the 
total number of swaps u^- for e^ and e^ for Uj , respectively, before reaching the 
solution. The Iterations column indicates the number of iterations (as defined 
in Section |4| the algorithm performed before reaching a solution. The Iter- 
ations with increases column shows the percentage of iterations where the 
number of changes increased from the previous iteration. We noticed that in the 
majority of iterations the number of changes decreased, which led to the quick 
convergence of the algorithm in the vast majority of cases. In all examples the 
starting point for the search was the ei, . . . , e„ base. The final column shows 
the percentage of problems where the algorithm was aborted due to going in a 
loop by allocating in some iteration a set of e^s and u^s that were encountered 
in a previous iteration. The percentage of loops was exponentially decreasing 
as the size increased and we did not observe any loops in any experiments on 
problem sizes of 30 or above. Overall, loops were observed in 0.1% of the exper- 
iments, so the heuristic algorithm was successful 99.9% of the time. A solution 
that we see for solving the problems that lead to a loop is to restart the algo- 
rithm from a different initial set of e^s and UjS. The problems ending in a loop 
were excluded from the detailed analysis that follows. 

^ htt p : // www .scilab. org/ 

^Any difference less than ±0.5 for integers and ±0.1 for percentages, respectively, is not 
shown, as deemed irrelevant for the analysis. 
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Table 1: Number of changes, iterations, iterations where the number of changes 
increased and loops for the various cone dimensions 



Size 


Changes 


Iterations 


Iterations with 


Loops 








increases [%] 


[%] 


2 


1 


1 





4.382 


3 


2 


1 


3.9 ±0.1 


4.278 


5 


4 


2 


13 ±0.2 


1.396 


10 


11 


3 


26 ±0.3 


0.273 


15 


17 


4 


30.3 ±0.3 


0.029 


20 


24 


4 


31.6 ±0.3 


0.007 


25 


30 


4 


31.2 ±0.3 


0.003 


30 


37 


4 


29.9 ±0.3 




50 


64 


5 


26.8 ±0.3 




75 


97 


5 


24.2 ±0.3 




100 


131 


5 


23.8 ±0.3 




200 


267 ± 1 


6 


19.9 ±0.8 




300 


409 ± 1 


6 


26.2 ±0.9 




500 


700 ±5 


7 


25.7±2.7 





More detailed analysis of the three main performance indicators of changes, 
iterations and iterations with increases was performed using boxplots as shown 
in Figures [T] [2] and [3j Although the total number of changes performed increases 
linearly with problem size (at a rate of less than 2 x n, even if considering 
maximum numer of changes) , this does not affect the performance substantially 
(see Figure [l] where the results were split into two parts for a clearer view) . 
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Figure 1: Boxplots of number of changes performed for the various cone dimen- 
sions 



The number of iterations is the crucial indicator of performance. As shown 
in Figure |2] the number of iterations reaches at most 11 (for sizes 15 and 20), 
but in 75% of cases has a value of 7 or below. We ran a few experiments on 
larger sizes, up to 1750, and the largest number of iterations we obeserved was 
13. Running experiments on very large problem sizes is problematic due to 
computer memory limitations and the Scilab built-in solving of linear systemsj^ 

^ Note that the time needed by one iteration substantially increases with problem size n 
as one iteration involves solving a linear system with n equations and n variables. 
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The major benefit of this heuristic algorithm is the small number of iterations 
even for very large number of cone dimensions. 
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Figure 2: Boxplots of number of iterations needed for the various cone dimen- 
sions 

As our heuristic algorithm seems to converge quickly, we wanted to know how 
frequently it deviates from the optimal path. An optimal path would consist of 
iterations with decreasing number of changes. Figure |3] shows the boxplots for 
the number of iterations where an increase in the number of changes took place. 
The maximum number of such iterations over all experiments is 4, but 75% of 
examples involved only one or no such iteration. This provides an explanation 
for the fast convergence: the algorithm very rarely deviates from the optimal 
path. 



6 Conclusion 

We presented a heuristic method of projection on simplicial cones based on 
Moreau's decomposition theorem. The heuristic algorithm presented in this 
note iteratively finds the projection onto a simplicial cone in a surprisingly 
small number of steps even for large cone dimensions in 99.9.% of the cases. 
We attribute the success to the fact that the algorithm rarely deviates from the 
optimal path, in every iteration it usually has to change less base values than 
in the previous iteration. We are planning to further extend the algorithm with 
random restart hoping to achieve 100% success rate. 
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Figure 3: Boxplots of number of iterations with increases in number of changes 
needed for the various cone dimensions 
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